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Abstract 

This  document  represents  the  final  report  on  the  various  scientific  activities  and  ac¬ 
complishments  relating  to  Grant  No.  FA9550-1 1-1-0194  over  its  period  of  performance, 

July  14,  2011  -  January  14,  2015.  The  project  had  the  following  five  overarching  tech¬ 
nical  objectives,  which  have  been  re-organized  under  three  broader  categories,  namely 
theoretical,  experimental,  and  computational  objectives  and  accomplishments  therein: 

1.  To  develop  physical  models  that  predict  how  surface  roughness,  texture,  and 
shape  alter  the  state  of  polarization  of  solar  illumination  scattered/reflected  by  a 
space  object; 

2.  To  develop  enhanced  but  approximate,  task-based  theoretical  performance-assessement 
tools  based  on  statistical  information  and  Bayesian  error  analysis  that  are  numer¬ 
ically  efficient,  polynomially  scalable,  and  generally  applicable  to  image  analysis, 
recovery,  and  reconstruction  tasks;  and  to  compare  the  actual  performance  of 
algorithms  with  theoretical  bounds  predicted  by  the  above  tools  and  to  explore 
ways  of  improving  performance  by  identifying  information  bottlenecks  and  weak 
links. 
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3.  To  design,  fabricate,  and  test  a  polarimetrically  enhanced  version  of  the  existing 
SD-CASSI  instrument  with  simple  polarization  sensitive  elements  added  to  the 
beam  train  to  generate  spectrally  compressive  polarimetric  data; 

4.  To  develop  and  benchmark  new  algorithms  for  highly  resolved  spatial  segmenta¬ 
tion,  material  identification,  shape  determination,  surface  characterization,  and 
extraction  of  rigid-body  kinematics  from  polarimetric  data  acquired  with  a  com¬ 
pressive  spectral  imaging  system  like  CASSI; 

5.  To  perform  an  exhaustive  analysis  of  a  newly  proposed  technique  for  the  recovery 
of  3D  shape,  rigid-body  kinematics,  and  surface  characterization  from  a  time 
series  of  2D  spectrally  compressive  polarimetric  images; 

Over  its  forty-two-month  performance  period,  the  research  grant  has  supported,  in 
part,  the  production  of  25  publications  (15  refereed  and  10  conference),  2  completed 
PhD  dissertations  at  UNM  (Srikanth  Narravula  -  defended  in  May  2014;  Rakesh  Kumar 
-  defended  on  April  2,  2015),  1  partially  completed  PhD  dissertation  at  Duke  (Tsung- 
Han  Tsai  -  supervised  by  Prof.  D.  Brady),  15  professional  presentations  (10  invited, 
5  contributed),  and  1  US  patent  application.  Based  on  research  partially  supported 
by  the  grant,  2  additional  manuscripts  are  currently  in  preparation  for  submission  to 
peer-reviewed  journals. 


1  Summary  of  Overall  Project  Accomplishments 

The  following  is  a  list  of  the  most  important  technical  accomplishments  of  the  project, 
separated  according  to  its  theoretical,  experimental,  and  computational  thrusts: 

1.1  Theoretical  Accomplishments 

1.1.1  Analysis  of  the  dependence  of  the  spatial  pBRDF  distribution  on  the  3D 
shape 

We  developed  a  unihed  theoretical  framework  to  analyze  compressive  sensing,  sparse  repre¬ 
sentation,  and  reconstruction  of  space  objects  in  the  combined  spatial-spectral  domain.  The 
usual  hyperspectral  imaging  approach  gathers  data  over  a  full  range  of  2D  spatial  coordinates 
and  ID  spectral  channels  in  what  is  called  a  data  cube.  In  most  cases,  such  data  typically 
carry  relatively  small  amounts  of  information  in  a  highly  redundant  manner,  since  space 
objects  whether  man-made  (like  satellites)  or  naturally  produced  (like  asteroids  and  other 
space  debris)  are  typically  composed  of  a  much  smaller  number  of  primitives,  e.g.,  spatially 
and  materially  homogeneous  elements  and  geometrically  simple  parts,  than  the  number  of 
voxels  in  the  full  data  cube. 

Theoretical  work  proceeded  along  two  different  lines.  The  first  entailed  constructing  sim¬ 
ple  forward  models  for  the  SOI  problem.  This  construction  begins  with  a  low-dimensional 
parametric  model  of  space-object  surfaces  using  a  superquadric  parameterization.  Superquadrics 
are  generalizations  of  ellipsiodal,  hyperboloidal,  and  toroidal  surfaces  based  on  changing  the 
powers  of  the  Cartesian-to-spherical  coordinate  relationships  for  such  ordinary  3D  surfaces. 


2 


Figure  1:  A  satellite  mock-up  composed  of  five  simple  superquadrics 


For  example,  a  super-ellipsoid  may  be  parameterized  as  follows: 

X  =ai  cos^^  ri\  cosa;|^^sgn(cosc<;); 

y  =02  cos^^  r]\  sina;|^^sgn(sina;); 

2;  =03!  sin  r]  I  ^^sgn  (sin  77),  (1) 

where  the  five  parameters  oi,  02, 03,  ei,  62  are  chosen  to  be  all  positive.  The  first  three  de¬ 
termine  the  spatial  scale  of  the  three  orthogonal  dimensions  in  the  object  while  the  last 
two  the  geometrical  shape  of  the  solid.  An  example  of  a  satellite  mock-up  composed  of  five 
different  superellipsoids  is  shown  in  Fig.  1.  By  parameterizing  the  relatively  small  number 
of  simple  geoemtrical  shapes  that  constitute  a  3D  satellite  body  via  superquadrics,  we  can 
turn  a  complex  surface  characterization  problem  involving  thousands  of  surface  pixels  into 
one  involving  a  rather  small  number  of  parameters.  These  parameters  can  be  estimated 
from  a  kinematic  sequence  of  2D  images  of  such  a  3D  body,  as  we  have  shown  previously  for 
low- dimensional  but  non-superquadric  parameterizations. 

Once  such  a  3D  space-object  surface  model  is  constructed,  it  may  be  ascribed  material 
properties,  specifically  what  material  comprises  what  superquadric  component  of  such  a 
body.  The  surfaces  are  then  given  texture.  To  do  so,  we  have  modeled  the  surfaces  of 
the  individual  superquadrics  as  being  microscopically  rough  so  they  do  not  scatter  sunlight 
simply  via  specular  reflection.  We  have  used  a  realistic  roughness  model  for  surfaces  in 
3D  that  may  be  described  statistically  by  means  of  a  distribution  of  planar  micro-facets 
with  varying  surface-normal  orientations  relative  to  the  underlying  idealized  microscopically 
smooth  surface.  Such  surfaces,  since  they  are  typically  made  from  dielectric  materials,  reflect 
and  refract,  so  there  is  both  a  single-scattering  specular  component  and  a  multiply-scattered 
diffuse  Lambertian  component.  It  is  the  latter  component,  which  gives  an  opaque  surface 
its  spectral  character,  i.e.,  color,  in  solar  reflection.  Both  these  components  contribute  to 
the  spectral  power  density  but  only  the  specular  albedo  has  a  polarimetric  content  that 
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(a)  (b) 

Figure  2:  The  2D  spatial  distribution  of  difference-polarization  intensity  in  solar  reflection 
from  a  sphere  for  four  different  solar  angles  for  (a)  smooth  surface;  (b)  rough  surface. 


carries  information  about  both  the  local  shape  and  roughness  of  the  space  object.  Thus, 
when  supplemented  with  a  sensor  noise  model,  a  model  for  atmospheric  propagation  and  a 
blur  model  for  the  observing  instrument,  our  calculations  allow  us  to  simulate  the  spectral 
and  polarimetric  content  of  any  observations  one  might  make  of  the  solar  albedo  from  a 
space  object  of  varying  local  shape  and  roughness  profile.  It  is  this  end-to-end  forward 
model  that  we  are  now  beginning  to  simulate  to  demonstrate  the  validity  of  our  compressive 
spectral-polarimetric  sensing  and  reconstruction  approach. 

The  power  spectral  density  (PSD)  of  the  diffuse  solar  albedo  depends  only  on  the  disper¬ 
sion  of  the  materials  comprising  the  surface,  rather  than  its  shape  or  roughness,  and  has  been 
extensively  measured  for  a  host  of  materials  commonly  present  on  a  satellite.  We  therefore 
do  not  focus  on  the  diffuse  PSD  in  our  effort,  but  analyze  the  specular  albedo  which,  as 
we  have  noted,  depends  not  just  on  material  dispersion  but  also  on  the  surface  shape  and 
roughness.  We  present  some  of  our  results  of  this  analysis  in  Fig.  2  where  we  display  the 
difference-polarization  (H-V)  images  of  a  roughened  sphere  of  Si  at  a  specific  wavelength 
and  four  different  solar  angles  of  illumination.  It  is  clear  that  roughness,  as  measured  by 
the  ratio  of  the  standard  deviation  of  the  heights  of  the  microfacets  and  its  length  affects 
the  footprint  of  the  polarimetric  signature  of  the  surface.  For  highly  smooth  surfaces,  the 
polarimetric  difference  image  is  limited  to  the  glint  spot,  but  as  the  surface  becomes  rougher, 
the  spot  containing  polarimetric  information  spreads  out  in  a  way  that  is  a  measure  of  the 
roughness  (for  fixed  shape),  as  we  see  in  Fig.  2(b). 

In  Fig.  3,  we  plot  the  depedence  of  the  angle-integrated  power  in  the  specular  component 
as  a  function  of  the  solar  angle  and  varying  roughness.  As  roughness  increases,  the  specular 
component  decreases  in  its  fractional  power  when  compared  to  the  diffuse  albedo,  so  the 
surface  becomes  more  Lambertian  and  acquires  its  spectral  signature  that  gives  it  a  color 
and  spectral  content  quite  distinct  from  the  white- light  spectrum  of  the  illuminating  sunlight. 
By  parameterizing  the  relatively  small  number  of  simple  geoemtrical  shapes  and  materials 
that  constitute  the  space  object  and  their  roughness  levels,  represented  by  the  ratio  of  the 
statistical  scatter  of  microfacet  heights  to  their  lateral  dimensions,  we  thus  have  an  excellent 
forward  model  for  varying  local  shapes,  spectral  reflectances,  and  roughness  levels  of  a  general 
space  object. 
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Reflection  from  Microfaceted,  Rough  Silicon  Surface 


Figure  3:  The  total  specular  return  fraction,  sHDR,  in  solar  reflection  vs.  the  solar  angle  for 
different  roughness  levels 

1.1.2  Extensions  of  the  micro- facet-based  roughness  model  to  include  hemi¬ 
spherical  pits  and  proposal  of  a  more  general,  composite  roughness  model 

For  ideally  smooth,  reflecting  objects,  the  solar-reflected  signal  is  rather  uninteresting,  char¬ 
acterized  by  simple  glint  points  at  which  the  condition  for  specular  reflectance  is  met  pre¬ 
cisely.  Such  objects  have  furthermore  spectral  traces  that  closely  reproduce  the  solar  spec¬ 
trum  which  they  reflect  rather  faithfully,  and  thus  appear  white.  But  such  idealized  surfaces 
are  just  that,  and  in  reality  they  undergo  processes  of  weathering  by  solar  wind  and  other 
particulate  surface  damage,  e.g.,  from  microdebris  flying  around,  that  roughen  them.  Rough¬ 
ening  spreads  the  reflected  radiance  over  the  surface  around  the  idealized  fglint  points  in  a 
way  that  depends  crucially  on  the  detailed  surface  shape  in  the  vicinity  of  the  glint  points. 
Roughening  can  thus  enable  surface  shape  extraction  if  a  detailed  physical  model  can  be 
derived  that  can  predict  such  dependence  on  shape  and  roughness.  Derivation  of  such  a 
forward  model  has  been  a  major  emphasis  of  our  theoretical  work.  We  previously  devel¬ 
oped  a  micro-facet-based  roughening  model  [5]  in  which  the  the  strength  of  the  polarimetric 
BRDF  from  any  point  on  the  object  surface  depends  on  its  surface  texture  as  well  as  the 
local  shape  and  orientation  relative  to  the  observation  direction  via  the  spatial  distribution 
of  the  reflected  power  around  specularly  refecting  glint  points. 

We  have  recently  extended  the  roughening  model  to  include  a  statistical  distribution  of 
hemispherical  pits,  rather  than  -  or  even  in  addition  to  -  microfacets,  from  which  surface 
returns  may  be  computed  semi-analytically.  These  pits  cause  both  included  and  excluded 
reflections  since  the  curved  surface  of  the  pits  allows  solar  light  to  be  reflected  in  a  whole 
range  of  directions  not  available  to  a  smooth  surface  without  pits,  while  the  same  curvedness 
also  excludes  points,  those  that  he  in  the  shadow,  from  being  operative  in  reflections  along 
a  corresponding  range  of  directions.  Such  included  and  excluded  reflections  depend  on  the 
depth-to-radius  ratio  of  the  pit,  and  thus  are  statistically  random  for  random  but  static 


5 


(a)  (b) 

Figure  4:  The  2D  angular  histogram  of  reflected  rays,  including  its  specular  and  diffuse  com¬ 
ponents,  for  solar  reflection  from  a  superquadric  for  two  different  solar  incidence  angles,  30° 
and  85°,  for  a  “pitted”  rough  surface  with  depth-to- radius  ratio  ranging  uniformly  between 
0  and  1. 


fluctuations  of  this  ratio  in  a  macroscopically  small  region  of  the  surface.  The  nature  of  the 
reflected  radiance,  as  captured  by  the  distribution  of  allowed  reflected  ray  directions,  thius 
changes  dramatically  when  compared  to  the  micro-facet-based  roughness  model.  This  is 
shown  in  the  following  figure  where  the  emergence  of  a  two-peaked  structure  in  the  reflected- 
ray  distribution  -  and  thus  in  the  BRDF  -  from  a  pit-based  model  is  easily  see,  in  sharp 
contrast  with  the  previous  model. 

We  can  also  separate  out  the  specular,  polarized  part  from  the  total  polarimetric  BRDF, 
so  we  can  determine  what  fraction  of  the  illumination  power  is  returned  in  the  unpolarized 
diffuse,  Lambertian  component.  This  requires  determining  the  singly-reflected  rays  of  all 
reflected  rays,  the  former  being  those  that  leave  the  pit  without  encountering  it  again.  This  is 
a  simple  geometrical  calculation  we  have  recently  completed.  We  are  currently  implementing 
this  pit-based  roughening  model  on  curved  surfaces  to  determine  how  for  such  surfaces  the 
polarimetric  BRDF,  like  in  the  micro-facet  model,  can  potentially  encode  the  3D  surface 
shape  and  surface  texture. 

We  are  using  our  two  models  of  roughening  to  construct  a  composite  roughness  model 
that  includes  these  two  different  roughness  primitives  -  microfacets  and  hemispherical  pits  - 
occurring  for  a  surface  exposed  to  the  two  different  mechanisms  of  roughening.  The  relative 
weight  of  the  two  BRDFs  corresponding  to  the  two  mechanisms  that  contribute  to  the  overall 
observable  BRDF  is  a  fitting  parameter  that  can  be  extracted  by  solving  the  inverse  problem 
of  reconstruction  and  parameter  recovery  based  on  our  forward  roughening  models.  Indeed, 
this  composite  roughness  model  is  easily  extended  further  to  include  a  whole  variety  of 
surface  roughening  and  damage  mechanisms  with  a  set  of  probability  weights  to  simulate 
realistically  roughened  surfaces  in  the  space  environment. 
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1.1.3  Developed  a  rotating  PSF  concept  [6,  7,  8]  to  perform  3D  imaging,  being 

extended  currently  to  include  polarimetric  imaging,  in  a  single  snapshot 

The  PI  proposed  and  developed  [6,  7,  8]  a  simple  approach  based  on  the  use  of  a  properly 
designed  pupil  phase  profile  that  is  divided  into  Fresnel-type  zones,  each  with  an  integral 
number  of  complete  phase  windings  that  changes  from  one  zone  to  the  next  by  a  fixed  integer, 
to  create  a  3D  point-spread  function  (PSF)  that  rotates  with  changing  defocus,  while  keeping 
its  transverse  shape  approximately  invariant.  This  shape  and  size  invariance  of  the  rotating 
PSF  can  easily  extend  over  ±  3-4  waves  of  defocus,  with  fairly  compact  dimensions  that 
can  provide  excellent  3D  precision  for  source  localization  and  imaging  of  3D  shapes.  Unlike 
Gauss-Laguerre  mode-based  approaches,  it  provides  highly  extended  depth  of  field  over  a 
large  transverse  field  of  view  to  enable  acquisition  of  large  3D  volumes.  The  off-center  spatial 
character  of  the  rotating  PSF  allows  it  to  resolve  with  ease  sources  that  are  in  the  same  line 
of  sight,  which  many  other  3D  imagers  cannot  typically  do  well.  Further,  it  generalizes 
readily  for  encoding  spherical  aberration  too  via  PSF  rotation. 

It  can  also  be  employed  to  effect  2D  digital  super-resolution  (DSR)  that  overcomes  the 
low- resolution  iFOV  limits  of  sensor  data  acquisition  by  collecting  multiple  image  data  frames 
that  are  different  only  in  the  axial  defocus  of  their  location.  Current  work  in  the  Pi’s  group  is 
beginning  to  address  the  prospects  of  DSR  by  replacing  the  more  traditional  sub-pixel-shifts 
by  changing  axial  defocus  of  PSF-rotated  image  frames. 

The  rotating-PSF  approach  can  be  extended  by  combining  its  spiral  phase  mask  with 
a  different  kind  of  phase  mask,  called  a  q-plate,  which  is  a  birefringent  phase  plate  with 
spiraling  optical  axis  of  phase  retardation.  Since  the  q-plate  can  convert  right  and  left 
circular  optical  polarizations  into  ±g  waves  of  orbital  angular  momentum,  this  combination 
mask  enables  one  to  generate  different  kinds  of  rotating  PSFs  for  different  source  polarization 
states.  In  other  words,  we  have  a  new  approach  to  perform  full  3D  polarimetric  imaging  in  a 
single  snapshot.  This  promising  concept,  with  tremendous  implications  for  SSA  and  general 
space-based  imaging  and  debris  localization,  has  been  proposed  for  further  study  under  a 
currently  pending  proposal  to  AFOSR  in  the  Imaging  Science  program  headed  by  Dr.  Julie 
Moses. 

1.1.4  Asymptotics  of  Bayesian  Error  Probability  and  Source  Super-Localization 

[9] 

We  carried  out  an  asymptotic  analysis  of  the  minimum  probability  of  error  (MPE)  in  inferring 
the  correct  hypothesis  in  a  Bayesian  multi-hypothesis  testing  (MHT)  formalism  using  many 
pixels  of  data  that  are  corrupted  by  signal  dependent  shot  noise,  sensor  read  noise,  and 
background  illumination.  We  performed  our  analysis  for  a  variety  of  combined  noise  and 
background  statistics,  including  a  pseudo-Gaussian  distribution  that  can  be  employed  to 
treat  approximately  the  photon-counting  statistics  of  signal  and  background  as  well  as  purely 
Gaussian  sensor  read-out  noise  and  more  general,  exponentially  peaked  distributions.  We 
subsequently  evaluated  both  the  exact  and  asymptotic  MPE  expressions  for  the  problem  of 
three-dimensional  (3D)  point  source  localization.  We  focused  specifically  on  our  recently 
proposed  rotating-PSE  imager  and  compared,  using  the  MPE  metric,  its  3D  localization 
performance  with  that  of  conventional  and  astigmatic  imagers  in  the  presence  of  background 
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and  sensor-noise  fluctuations. 

Based  on  the  powerful  asymptotic  error  analysis,  we  determined  that  the  rotating  PSF 
imager  has  the  best  super-localization  performance,  of  all  3D  imagers  considered  here,  over 
a  large  focal  volume  containing  multiple  point  sources.  This  work  was  also  presented  at 
the  2014  International  Quantitative  Bio-Imaging  workshop  as  a  poster,  and  elicited  much 
interest  from  the  single-molecule  bio-imaging  community. 

The  analysis  represents  the  achievement  of  an  important  theoretical  milestone  in  error 
analysis  of  single-molecule  localization  that  is  not  limited  in  applicability  as  are  those  based 
on  the  Fisher  information  and  associated  Cramer-Rao  bound  on  the  error  of  estimation. 

1.1.5  Asymptotics  of  Bayesian  Error  Probability  and  2D  Pair  Superresolution 

[10] 

We  applied  the  asymptotic  Bayesian  multi- hypothesis  testing  (MHT)  based  error  analysis 
developed  in  the  previous  paper  (see  above)  to  treat  the  problem  of  superresolution  imaging 
of  a  pair  of  closely  spaced,  equally  bright  point  sources.  The  analysis  exploits  the  notion 
of  the  minimum  probability  of  error  (MPE)  in  discriminating  between  two  competing  equi- 
probable  hypotheses,  a  single  point  source  of  a  certain  brightness  at  the  origin  vs.  a  pair 
of  point  sources,  each  of  half  the  brightness  of  the  single  source  and  located  symmetrically 
about  the  origin,  as  the  distance  between  the  source  pair  is  changed.  For  a  Gaussian  point- 
spread  function  (PSF),  the  analysis  makes  predictions  on  the  scaling  of  the  minimum  source 
strength,  expressed  in  units  of  photon  number,  required  to  disambiguate  the  pair  as  a  function 
of  their  separation  in  both  the  signal-dominated  and  background-dominated  regimes.  Certain 
logarithmic  corrections  to  the  quartic  scaling  of  the  minimum  source  strength  with  respect 
to  the  degree  of  superresolution  characterize  the  signal-dominated  regime,  while  the  scaling 
is  purely  quadratic  in  the  background-dominated  regime.  For  the  Gaussian  PSF,  general 
results  for  arbitrary  strengths  of  the  signal,  background,  and  sensor  noise  levels  were  also 
presented. 

The  applicability  of  this  error  analysis  to  more  complicated  imagers  is  being  currently 
pursued.  Our  hope  is  to  provide  a  comprehensive  Bayesian  inference  based  description  and 
quantitative  understanding  of  other  imagers  in  performing  3D  superresolution  imaging  of 
closely  spaced  point  sources. 

1.1.6  Bayesian  Error  Based  Sequences  of  Statistical  Information  Bounds  [11] 

The  relation  between  statistical  information  and  Bayesian  error  has  been  sharpened  by  de¬ 
riving  finite  sequences  of  upper  and  lower  bounds  on  the  equivocation  entropy  (EE)  in  terms 
of  the  minimum  probability  of  error  (MPE)  and  related  Bayesian  quantities.  The  well  known 
Eano  upper  bound  and  Eeder-Merhav  lower  bound  on  the  EE  are  tightened  by  including  a 
succession  of  posterior  probabilities  starting  at  the  largest,  which  directly  controls  the  MPE, 
and  proceeding  to  successively  lower  ones.  A  number  of  other  interesting  results  were  also 
derived,  including  a  sequence  of  upper  bounds  on  the  MPE  in  terms  of  a  previously  intro¬ 
duced  sequence  of  generalized  posterior  distributions.  The  tightness  of  the  various  bounds 
was  numerically  evaluated  for  a  simple  example.  These  new  bounds  are  not  just  compu¬ 
tationally  efficient  but  they  are  sufficiently  simple  to  interpret  analytically,  and  can  thus 


provide  excellent  framework  in  which  to  analyze  information  processing  systems  that  use 
finite  codewords  for  which  some  of  the  more  powerful  asymptotic  approaches  cannot  be  used 
without  incurring  large  errors. 

Our  plan  for  current  and  future  work  is  to  employ  some  of  these  bounds  to  analyze 
fundamental  upper  limits  of  performance  of  3D  imaging  systems  like  rotating-PSF  imagers 
when  only  a  few  pixels  of  image  data  are  available  or  practical,  e.g.,  in  the  point-source 
localization  scenario  where  many  sources  cover  a  relatively  small  image  area  so  only  a  few 
pixels  per  source  can  be  used  to  localize  each  source.  An  asymptotic  analysis  would  not  be 
sufficiently  accurate  in  such  a  situation,  but  the  bounds  developed  in  this  work  can  be  used 
with  much  accuracy. 

1.2  Experimental  Accomplishments 

1.2.1  Coded  aperture  snapshot  spectral  polarization  imaging  -  early  version 
using  two  birefringent  crystals  in  combination  with  a  single  disperser  and 
a  grayscale  coded  mask 

We  demonstrated  a  coded  aperture  snapshot  spectral  polarization  imager  based  on  com¬ 
pressive  sensing.  In  a  snapshot,  this  system  can  calibrate  the  difference  of  the  two  linear 
orthogonal  polarization  states  along  with  their  spectral-spatial  signature  with  large  field  of 
view.  A  grayscale  2D  detector  array  is  utilized  to  record  the  coded  mixture  of  spatial,  spec¬ 
tral  and  polarization  information  of  the  object.  In  the  end,  a  fast  imaging  reconstruction 
algorithm  called  TwIST  is  utilized  to  reconstruct  the  polarized  spectral  image. 


Figure  5:  System  setup. 

Figure  5  is  a  photograph  of  the  experimental  prototype.  The  objective  lens  in  this  setup 
is  a  commercialized  objective  form  Jenoptik  (Easthampton,  Massachusetts).  The  coded 
aperture  is  a  lithography  patterned  on  a  quartz  substrate  with  an  anti-reflection  coating  on 
both  sides,  with  1988*1988  elements  random  binary  pattern  that  are  each  7.4  jim  square. 
The  birefringent  material  are  two  calcite  crystals  which  manufactured  by  United  Crystals 
(Port  Washington,  New  York),  with  8mm  thickness,  2.5  cm  square  clear  aperture  and  45° 
optical  axis  angle  relative  to  the  light  propagation  axis.  The  magnification  of  the  relay 
lens  is  0.38.  Finally,  the  detector  is  an  8-bit  monochrome  CMOS  camera  form  Aptina 
(San  Jose,  California)  with  4384*3288  pixels  that  are  each  1.4  jim  square.  All  of  these 
system  components  are  aligned  on  a  mounting  plate.  The  detector  and  the  objective  lens  are 
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carried  by  translation  stages  in  order  to  provide  the  ability  to  adjust  the  focusing  distance. 
The  detector  was  placed  on  a  rotation  mount  to  ensure  the  proper  alignment  of  the  beams 
dispersion  direction  across  the  image  plane. 

Figure  6  represents  the  design  and  the  dispersion  relation  of  the  birefringent  crystal  pair. 
The  optical  axes  of  the  first  and  the  second  calcite  plates  orientate  45°  relative  to  the  y-z 
plane  and  x-z  plane,  respectively.  Therefore,  the  ordinary  ray  (o-ray)  of  the  first  birefringent 
plate  will  become  the  extraordinary  ray  (e-ray)  of  the  second  plate  and  be  shifted  parallel  to 
its  original  propagation  direction.  This  image  displacing  effect  split  two  orthogonal  polarized 
images.  In  addition,  the  displacing  distance  follows  the  refraction  law,  which  results  in  image 
sharing  at  wavelength  dependent  location.  When  the  images  pass  through  the  crystal  pair, 
those  birefringent  crystals  can  generate  the  polarization  selective  image  displacement  and 
wavelength  dependent  image  dispersion.  Both  modulations  are  due  to  the  double  refraction 
effect  of  the  calcite  crystals. 


Figure  6:  Design  of  the  birefringent  crystals. 

On  the  detector  plane,  the  measurement  in  this  system  is  a  two-dimensional  projection 
of  two  separated  three-dimensional  spatial-spectral  data  cube  which  include  the  information 
form  the  scene.  The  data  reconstruction  form  this  compressive  sensing  based  on  assuming 
the  object  has  the  sparsity  in  some  basis,  which  making  the  object  highly  compressible. 
Here  the  reconstructing  process  relies  on  the  two-step  iterative  shrinkage  and  thresholding 
(TwIST)  algorithm.  We  also  assume  the  piecewise  smoothness  of  object  in  spatial  domain 
and  the  sparsity  in  the  wavelet  basis  is  enforced  by  utilizing  the  total  variation  regularization. 

To  test  the  performance  of  the  system  and  the  reconstruction  algorithm,  we  measured 
negative  1951  USAF  resolution  test  chart  to  validate  the  ability  of  the  TwIST  algorithm 
to  estimate  the  variation  in  spectral  signature  and  partial  polarization  state  distribution  in 
an  image.  The  test  chart  was  illuminated  by  a  Tungsten  light  bulb  and  filtered  by  a  linear 
polarizer.  The  fight  was  guided  by  a  multimode  fiber  and  collimated  by  a  commercialize 
beam  expander.  A  linear  polarizer  was  placed  after  the  collimator  to  control  the  incident 
polarization  state  of  the  system.  The  transmission  image  of  the  test  chart  was  modulated  and 
then  acquired  by  the  measurement  system.  Figure  7  shows  a  detector  measurement  of  the  test 
chart.  The  dispersion  and  splittering  which  was  generated  by  the  two  birefringent  crystals 
redistributed  the  intensity  of  the  incident  image.  As  shown  in  the  figure,  the  measurement 
on  the  detector  plane  is  the  superposition  of  the  two  overlapping  images;  each  image  is 
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Figure  7:  Detector  measurement. 


modulated  by  the  coded  pattern.  A  400-680  nm  band-pass  filter  was  mounted  in  front  of  the 
objective  lens  to  confine  the  measurement  of  the  scene  within  the  corresponding  wavelength 
channels. 

Figure  8  presents  the  reconstruction  of  19  wavelength  channel  and  two  Stokes  vectors 
between  400  and  680  nm.  The  left  part  of  the  reconstruction  shows  the  So  term  of  the 
scene,  which  is  total  irradiance  of  the  object;  the  right  part  shows  the  6*1  term,  which 
represents  the  polarization  difference  in  intensity.  The  spatial  modulation  pattern  and  the 
image  overlapping  in  Figure  5  had  been  removed  and  clear  resolution  charts  with  sharp  edge 
had  been  reconstructed  by  the  TwIST  algorithm  in  majority  of  the  wavelength  channels. 
The  brightness  in  the  reconstructed  figure  represents  the  relatively  intensity  between  the 
wavelength  channels.  Since  the  azimuth  angle  of  the  incident  polarizer  was  30°  relative  to 
the  plane  of  the  optical  table,  the  ideal  value  of  the  Si  terms  is  0.5.  Based  on  the  result 
in  Figure  4,  the  brightness  of  the  images  in  all  channels  are  around  half  of  value  as  the 
corresponding  irradiance  channels,  which  fits  the  theoretical  value. 

To  test  the  performance  of  the  spectral  reconstruction,  we  compared  the  reconstructed 
spectral  signature  of  the  test  chart  and  the  reference  spectrum  was  acquired  using  a  com¬ 
mercial  Ocean  Optics  USB2000  spectrometer.  For  comparison,  the  continuous  reference 
spectrum  was  integrated  into  the  19  spectral  bands  based  on  the  channel  width  of  this 
system.  Figure  9  shows  the  reconstructed  spectrum  of  the  test  chart  in  blue  star  and  the  ref¬ 
erence  spectrum  in  red  line.  Both  spectrums  were  normalized  to  the  maximum  value  in  their 
curves.  The  reconstructed  spectrum  represents  the  spectral  distribution  of  the  average  in¬ 
tensity  within  a  selected  area  in  the  test  chart.  The  comparison  demonstrates  the  agreement 
between  the  reconstructed  and  the  reference  spectrum,  indicating  that  the  reconstruction 
can  correctly  identify  the  spectral  data. 
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Figure  8:  Reconstruction  result  of  a  negative  1951  USAF  resolution  test  chart. 


1.2.2  Coded  aperture  snapshot  spectral  polarization  imaging  -  improved  version 
using  coded  aperture,  polarizers,  and  LC-SLM  to  create  joint  polarization 
and  spectral  encoding 

We  next  employed  a  vastly  improved  design  based  on  the  use  of  a  liquid  crystal  on  silicon 
modulator  to  jointly  code  spatial,  spectral  and  polarization  features  for  snapshot  spectral 
imaging.  We  demonstrate  compressive  sampling  of  megapixel  multispectral  image  on  a  two- 
dimensional  (2D)  detector  array.  The  reconstruction  recovers  the  encoded  2D  measurement 
to  a  4D  data  cube  The  4D  data  cube  includes  megapixels  in  spatial  domain,  15  wavelength 
channels  between  450  and  680  nm  and  their  horizontal  and  vertical  polarized  irradiance. 

A  schematic  of  this  SLM  based  CASSI  system  is  shown  in  Fig.  10.  An  objective  lens  LI 
images  the  scene  from  a  remote  distance,  and  the  lenses  L2,  L3  relay  the  scene  through  the 
beam  splitter  and  the  wave  plate  onto  the  modulator.  The  SLM  modulates  the  scene  with 
spatially  distributed  polarization  state  rotation  and  then  reflects  the  image  back.  The  SLM 
generated  degree  of  polarization  state  rotation  is  a  wavelength  dependent  value.  The  imaging 
lenses  L3  and  L4  relay  the  modulated  image  onto  the  detector  plane  through  the  polarizer. 
To  achieve  best  contrast  in  every  wavelength,  we  insert  an  achromatic  quarter  wave-plate 
in  the  light  path  in  order  to  compensate  the  extra  phase  retardation  in  the  liquid  crystal 
device.  Finally,  the  detector  array  records  modulated  patterns  in  every  wavelength,  which 
involves  the  spectral  and  spectral,  and  polarization  modulated  information  of  the  object. 

Fig.  11  is  a  photograph  of  the  experimental  setup.  The  setup  includes  a  60  mm  objective 
lens  (Jenoptik),  one  70  mm  achromatic  lens  (Newport),  a  dichroic  beam  splitter  (Newport), 
an  achromatic  quarter  wave  plate,  a  liquid  crystal  based  SLM  (Pluto,  Holoeye),  two  25  mm 
imaging  lenses  (Computar),  a  sheet  polarizer (CVI)  ,and  a  monochrome  CCD  camera  (Pike 
F-421,  AVT)  with  2048  x  2048  pixels  that  are  7.4  /rm  square.  A  band  pass  Alter  (Badder) 
is  mounted  on  the  objective  lens  to  cut  the  ultra-violet  and  infrared.  The  transmission  of 
the  filter,  the  reflection  of  the  SLM  and  the  quantum  efficiency  of  the  detector  limit  the 
wavelengths  of  interest  of  this  system  between  450  nm  and  680  nm.  The  scattered  image  of 
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Figure  9:  Comparison  between  the  reconstructed  and  the  reference  spectrum. 


the  scene  is  first  collected  by  the  objective  lens  then  projected  on  the  SLM.  After  modulated 
by  the  SLM  with  transmission  patterns,  the  scene  is  reflected  and  then  propagated  to  the 
color  insensitive  detector. 

The  SLM  plays  an  important  role  in  this  system  because  this  revised  design  of  polar¬ 
ization  spectral  imager  replaces  both  the  coding  device  and  the  dispersive  element  by  this 
single  modulator.  This  modulator  is  based  on  a  reflective  Liquid  Crystal  on  Silicon  (LCoS) 
microdisplay  technique,  which  has  a  1920x1080  pixels  active  area  with  8  /rm  pixel  pitch. 
Being  down  magnified  by  two  imaging  lenses,  the  width  of  each  pixel  on  the  SLM  becomes 
7.4  jim  square.  Considering  a  perfect  alignment,  it  is  possible  to  map  each  detector  pixel 
with  one  corresponding  pixels  on  the  SLM.  Each  pixel  on  the  microdisplay  can  be  indepen¬ 
dently  addressed  by  an  8-bit  voltage  value  to  change  the  amount  of  phase  retardation  to  the 
incident  polarized  light.  After  propagating  back  through  the  polarizer  splitter,  the  amount 
of  retardation  for  two  orthogonal  polarized  incident  becomes  different  transmission  to  every 
voxel  in  the  image  data  cube  corresponding  to  its  spatial  voltage  address  and  wavelength. 
This  characteristic  will  be  utilized  to  encode  the  image  in  spatial  and  spectral  domains. 

Here  we  used  a  random  8-bit  (0-255)  gray  scale  voltage  code  on  the  SLM  to  generate 
a  group  of  different  transmission  patterns  for  each  wavelength  slice  (450  to  680  nm)  in  the 
image  data  cube.  The  voltage  address  varies  from  0  to  255,  which  represents  0  to  3  tt  phase 
delay  at  630  nm.  The  transmission  to  voltage  address  and  wavelength  can  be  estimated 
by  calibrating  the  SLM.  The  calibration  process  records  the  transmission  under  the  two 
orthogonally  polarized  monochromatic  illuminations  and  homogeneously  applied  voltage  on 
the  modulator.  It  also  includes  the  transmission  of  all  the  optics,  the  quantum  efficiency  of 
the  detector,  and  the  phase  modulation  of  the  SLM.  We  averaged  the  transmission  across  the 
center  area  of  the  SLM  and  calculate  its  amplitude  modulation,  which  is  shown  in  Fig.  12. 
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Figure  10:  The  system  schematic.  LI  is  the  objective  lens  images  the  scene.  L2  and  L3 
relay  the  image  onto  the  SLM.  L3  and  L4  project  the  encoded  image  on  the  detector.  The 
polarizer  is  used  to  achieve  intensity  modulation  mode  of  the  SLM.  The  quarter  wave  plate 
compensates  phase  retardation  of  the  SLM  to  achieve  best  contrast. 


Based  on  the  calibration,  within  450  nm  and  680  nm,  an  8-bit  voltage  variation  on  the  SLM 
generates  a  gray  scale  transmission  ratio  between  0.08  and  0.96. 

Here  the  liquid  crystal  SLM  generates  joint  wavelength  and  polarization  dependent  mod¬ 
ulation  for  each  spectral  plane  of  the  scene.  The  modulated  spectral  slices  are  integrated  by 
a  2D  detector,  but  the  spectral  information  in  the  4D  data  cube  can  be  distinguished  based 
on  the  corresponding  transmission  patterns  between  all  wavelength  channels. 

1.2.3  A  Further  Improved  Version  of  the  SLM-CASSPI  Instrument  with  Full 
Linear  Polarimetric  Coding 

We  decompose  the  compressed  signal  into  0*^,  90*^,  45°,  and  135°  polarization  channels  for  lin¬ 
ear  polarization  state  estimation.  This  linear  polarimetry  would  satisfy  several  applications 
without  significant  circular  polarizations,  such  as  most  of  the  natural  scenes. 

The  LC  cell  has  a  property  of  electrical  controllable  optical  anisotropy,  which  enables 
the  signal  modulation  for  this  camera.  The  long  axis  and  the  short  axis  of  each  layer  of 
the  LC  cell  provide  different  refractive  indices  to  the  incident  wave.  By  controlling  the 
orientation  of  the  LC  cell  through  tuning  the  applied  voltage,  the  SLM  can  be  considered 
as  a  variable  wave-plate.  We  sampled  the  vertical  fraction  of  the  irradiance  to  transfer  the 
phase  modulation,  given  by  jS,  into  a  detector  recognizable  amplitude  modulation,  which  can 
be  described  as 

/  =  (1/2)  [(5*0  —  >5i  cos (2/3)  S 2  sin (2/3)] , 

in  terms  of  the  first  three  Stokes  parameters,  Sq,  Si,  S2. 

Fig.  14  shows  the  amplitude  modulation  provided  by  the  different  polarization  and  color 
channels.  Each  sub-figure  describes  the  relationship  between  the  transmission  pixel  counts 
and  the  applied  voltage  on  the  SLM.  We  note  that  this  camera  has  freedom  of  choosing 
the  basis  of  polarization  state  decomposition.  Here  we  use  linear  horizontal,  linear  vertical, 
linear  45,  and  linear  135  as  the  decomposing  basis  to  analyze  all  linear  polarization  states. 

Fig.  11  shows  the  experimental  prototype  of  the  compressive  camera.  A  series  of  lenses 
projects  the  scene  on  the  SLM,  which  is  the  intermediate  image  plane.  A  random  voltage 
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Figure  11:  The  experimental  prototype  of  the  system. 


pattern  is  applied  on  the  SLM,  which  provides  a  spectral  and  polarization  dependent  ampli¬ 
tude  modulation  to  encode  the  object.  After  reflected  by  the  silicon  layer  of  the  SLM  and 
then  Altered  by  the  polarizer,  the  applied  phase  retardation  becomes  a  wavelength  and  po¬ 
larization  dependent  amplitude  modulation  which  is  added  to  the  scene.  We  apply  a  quarter 
wave-plate  as  a  compensator  to  increase  the  contrast  ratio  of  the  coding.  The  modulated 
image  will  then  recorded  by  the  color  detector.  Based  on  the  design,  this  camera  has  a  25 
field  of  view  and  a  12.5  mm  clear  aperture.  The  detector  has  a  1280x720  spatial  resolution 
with  4.08um  square  pixel  size.  The  modulation  process  can  be  represented  by  the  following 
mathematical  model: 


9mn 


[/||  (x,  y,  A)7j|  (x,  y,  X)+f±ix,  y,  A)rj_(x,  y,  A)]  rect  (x/ A-m,  y / A-n)dxdydX+Wr, 


where  g  represents  the  power  spectral  density  in  the  detector  plane,  /  the  spatial-spectral 
distribution  of  the  scene,  and  T±  and  Tj|  are  the  wavelength  and  polarization  dependent 
transmission  code  patterns.  We  note  that  the  size  of  the  detector  pixel  is  A  and  its  corre¬ 
sponding  pixel  noise  is  denoted  by  w.  We  depict  the  4D  discrete  sampling  function  in  matrix 
form  as  g  =  Hf  -|-  w,  where  H  represents  the  forward  matrix  of  the  system,  f  the  vectorized 
and  discretized  object  datacube,  and  w  represents  the  sensor  noise  vector.  The  H  matrix 
approximates  the  sampling,  encodes  the  color  and  polarization  channels  to  compress  the  4D 
color  and  polarization  datastream  f  on  measurement  g.  We  experimentally  calibrate  the 
forward  matrix  H  by  illuminating  the  camera  under  four  identical  polarization  states.  We 
assign  linear  horizontal,  linear  vertical,  linear  45°,  and  linear  135°  polarization  channels  as 
the  identical  states.  The  light  source  we  use  is  a  white  light  LED  filtered  by  a  rotatable 
polarizer.  We  record  the  modulation  patterns  of  these  four  polarizations  sequentially.  This 
calibration  considers  the  errors  which  might  break  the  ideal  mapping  relationship,  such  as 
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Figure  12:  The  intensity  modulation  of  the  SLM  for  horizontal  polarized  light.  Given  dif¬ 
ferent  voltages  and  incident  wavelengths,  the  SLM  combined  with  the  optics  generates  cor¬ 
responding  transmission  code. 


the  aberration  of  the  optics  and  the  sub-pixel  misalignment. 

Figs.  15  and  16  show  the  experimental  result  of  the  compressive  camera.  The  objects 
are  polarization  filtered  color  bricks  and  a  scene  of  a  parking  lot.  Both  of  them  include  the 
spatial,  spectral  and  polarization  complexity.  Both  scenes  are  reconstructed  by  Two-step 
Iterative  Shrinkage  Thresholding  (TwIST)  algorithm  combined  with  total  variation  (TV) 
regularization  function.  The  reconstruction  shows  the  linear  horizontal,  linear  vertical,  linear 
45'^,  and/or  linear  135'^  polarization  channels  with  red,  green,  and  blue  color  channels. 

This  single-shot  compressive  color  polarization  imager  presents  an  integration  of  color 
and  polarization  compression.  Such  camera  uniquely  encodes  and  decomposes  the  scene  into 
color  polarization  images  through  the  phase  modulation  of  a  SLM  along  without  any  dis¬ 
persive  element.  This  technique  eliminates  mechanical  movements  that  hinder  conventional 
polarimetry.  The  experimental  results  show  clear  spatial  resolution  and  low  noise  reconstruc¬ 
tions  in  the  experiments.  Future  applications  will  use  the  polarization  sensitivity  to  analyze 
the  surface  information,  such  as  the  curvature  and  the  roughness. 

A  similar  setup  can  be  used  in  compressive  spectral  imaging.  Based  on  the  singular 
value  decomposition,  this  SLM  modulation  provides  six  measurement  modes  in  spectral 
multiplexing.  Therefore,  it  can  effectively  multiplex  a  spectrum  which  has  up  to  six  major 
spectral  bands/peaks  in  the  sensing  region  (400  to  680  nm)  without  losing  information.  This 
criterion  can  satisfy  most  of  the  spectral  distribution  of  the  natural  scenes,  which  have  smooth 
spectrum  and  only  require  two  to  three  modes  to  fully  represent  their  major  distribution  in 
the  spectral  domain.  The  minor  detail  can  be  recovered  by  the  spectral  constraint  of  the 
reconstruction  algorithm,  such  as  interpolating  the  relative  spectral  density  to  correct  the 
reconstruction. 

Here  we  use  an  SLM  to  encode  the  3D  spatial-spectral  information  on  a  2D  gray-scale 
detector  -  a  similar  CDMA  strategy  to  CASSI  for  snapshot  hyperspectral  imaging.  The  SLM 
provides  wavelength-dependent  transmission  patterns  to  multiplex  every  spectral  channel  for 
the  compressive  measurement.  The  hyperspectral  slices  can  be  separated  by  using  inversion 
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Figure  13:  The  intensity  modulation  of  the  SLM  for  vertical  polarized  light.  Given  different 
voltages  and  incident  wavelengths,  the  SLM  combined  with  the  optics  generates  correspond¬ 
ing  transmission  code. 

algorithms.  We  calibrated  the  spectral  channels  under  a  5  nm  spectral  bandwidth  by  using 
a  monochromator  (iHR320,  Horiba). 

Figs.  17  and  18  show  the  preliminary  reconstruction  results  of  this  compressive  hyper- 
spectral  camera.  The  objects  are  color  candies  which  are  illuminated  by  a  tungsten  light 
source. 

1.3  Computational  Accomplishments 

1.3.1  Randomized  Singular- Valued  Decomposition  (SVD)  Methods  in  Hyper- 
spectral  Imaging  [12] 

We  advanced  and  tested  a  randomized  singular  value  decomposition  (rSVD)  method  for 
the  purposes  of  lossless  compression,  reconstruction,  classification,  and  target  de-  tection 
with  hyperspectral  (HSI)  data.  Recent  work  in  low-rank  matrix  approximations  obtained 
from  random  projections  suggest  that  these  approximations  are  well-suited  for  randomized 
dimensionality  reduction.  Approximation  errors  for  the  rSVD  were  eval-  uated  on  HSI  and 
comparisons  were  made  to  deterministic  techniques  and  as  well  as  to  other  randomized  low- 
rank  matrix  approximation  methods  involving  compressive  prin-  cipal  component  analysis. 
Numerical  tests  on  real  HSI  data  suggested  that  the  method  is  promising,  and  is  particularly 
effective  for  HSI  data  interrogation. 

1.3.2  Deblurring  and  Sparse  Unmixing  For  Hyperspectral  Images  [14] 

We  developed  a  new  approach  to  joint  deblurring  and  sparse  unmixing  of  hyperspectral 
images  based  on  a  total  variation  regularization  method: 

min  i  \\HXA  -Y\\%  +  pn  ||,Y||, ,  +  ,<2TV(A') 


Vertical 
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Figure  14:  The  amplitude  modulation  tests  for  three  color  channels  and  four  polarization 
channels.  The  horizontal  axis  represents  the  8-bit  applied  voltage  address  on  the  SLM.  The 
vertical  axis  is  the  average  pixel  count. 

where  H  is  an  n-by-n  hyperspectral  data  acquisition  blurring  matrix,  and  fii  and  fi2  are  two 
positive  regularization  parameters  used  to  control  the  importance  of  the  sparsity  term  and 
the  total  variation  term,  respectively.  Here 

n  m 

i=l  j=l 

In  the  model  we  incorporated  blurring  operators  for  dealing  with  blurring  effects  accross  the 
spectral  bands,  and  in  particular,  blurring  operators  for  hyperspectral  imaging  whose  PSFs 
are  generally  system  dependent  and  result  from  axial  optical  aberrations  in  the  acquisition 
system.  An  alternating  direction  method  was  developed  to  solve  the  resulting  optimization 
problem  efficiently.  According  to  the  structure  of  the  total  variation  regularization  and 
sparse  unmixing  in  the  model,  the  convergence  of  the  alternating  direction  method  can  be 
guaranteed.  Experimental  results  were  reported  to  demonstrate  the  effectiveness  of  the  TV 
and  sparsity  model  and  the  efficiency  of  the  proposed  numerical  scheme,  and  the  method 
was  compared  to  the  recent  Sparse  Unmixing  via  variable  Splitting  Augmented  Lagrangian 
and  Total  Variation  (SUnSAL-TV)  method  by  lordache,  Bioucas-Dias,  and  Plaza. 

Here  we  assumed  that  the  PSF  leading  to  the  blurring  matrix  for  hyperspectral  data 
has  already  been  estimated.  As  discussed  above,  hyperspectral  image  degradation  is  often 
“system  dependent,”  and  the  point-spread  function  (PSF)  of  the  image  acquisition  system 
needs  to  be  identified.  The  system  PSF  relates  generally  to  axial  optical  aberrations,  while 
spatial  PSFs  relate  to  the  usual  image  blurring  problems  such  as  defocus,  motion,  imaging 
through  a  medium,  etc.  Axial  optical  aberrations  can  lead  to  a  significant  blurring  of  image 
intensities  in  certain  parts  of  the  spectral  range.  These  axial  optical  aberrations  arise  from 
the  index  of  refraction  variations  that  is  dependent  on  the  wavelength  of  incident  light.  In 
a  recent  paper  by  Spiclin,  et  ah,  the  authors  assumed  a  model  of  the  PSF  to  be  a  “linear 
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Figure  15:  The  measurement,  the  references,  and  the  reconstruction  of  toys  which  is  filtered 
by  two  orthogonal  sheet  polarizers.  The  upper  left  shows  the  azimuth  angle  of  the  polarizer. 
The  upper  middle  and  upper  right  show  the  irradiance  of  the  object  and  the  raw  measurement 
of  the  compressive  polarization  camera,  respectively.  The  second  row  shows  the  reference 
images  measured  by  a  camera  filtered  by  a  rotatable  polarizer.  The  third  row  shows  the 
pseudo-color  reconstruction. 


combination  of  Gaussian  functions”  across  the  different  spectral  bands,  therefore  the  overall 
PSF  blur  identification  process  is  reduced  to  finding  only  the  corresponding  scalar  weight 
for  each  Gaussian  function. 

1.3.3  Sparse  Nonnegative  Matrix  Underapproximation  and  its  Application  to 
Hyperspectral  Image  Analysis  [15] 

Dimensionality  reduction  techniques  such  as  principal  component  analysis  (PGA)  are  power¬ 
ful  tools  for  the  analysis  of  high-dimensional  data.  In  hyperspectral  image  analysis,  nonneg¬ 
ativity  of  the  data  can  be  taken  into  account,  leading  to  an  additive  linear  model  called  non¬ 
negative  matrix  factorization  (NMF),  which  improves  interpretability  of  the  decomposition. 
Recently,  another  technique  based  on  underapproximations  (NMU)  has  been  introduced, 
which  allows  the  extraction  of  features  in  a  recursive  way,  such  as  PGA,  but  preserving 
nonnegativity,  such  as  NMF.  Moreover,  in  some  situations,  NMU  is  able  to  detect  automat¬ 
ically  the  materials  present  in  the  scene  being  imaged.  However,  for  difficult  hyperspectral 
datasets,  NMU  can  mix  some  materials  together,  and  is  therefore  not  able  to  separate  all  of 
them  properly.  In  this  paper  we  introduce  sparse  NMU  by  adding  a  sparsity  constraint  on 
the  abundance  matrix  and  use  it  to  extract  materials  individually  in  a  more  efficient  way 
than  NMU.  This  was  experimentally  demonstrated  on  the  HYDIGE  images  of  the  San  Diego 
airport  and  the  Urban  dataset. 
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Figure  16:  The  measurement,  the  references,  and  the  reconstruction  of  a  scene  of  a  parking 
lot.  The  upper  left  and  upper  right  show  the  irradiance  of  the  object  and  the  raw  mea¬ 
surement  of  the  compressive  polarization  camera,  respectively.  The  second  row  shows  the 
reference  images  measured  by  a  camera  filtered  by  a  rotatable  polarizer.  The  third  row  shows 
the  pseudo-color  reconstruction.  Notice  that  the  reflections  on  the  windows  and  on  the  rear 
screen  of  cars  are  polarized. 

1.3.4  Randomized  Methods  in  Lossless  Compression  of  Hyperspectral  Data  [16] 

We  evaluated  recently  developed  randomized  matrix  decomposition  methods  for  fast  lossless 
compression  and  reconstruction  of  hyperspectral  imaging  (HSI)  data.  The  simple  random 
projection  methods  have  been  shown  to  be  effective  for  lossy  compression  without  severely 
affecting  the  performance  of  object  identification  and  classification.  We  built  upon  these 
methods  to  develop  a  new  double-random  projection  method  that  may  enable  security  in  data 
transmission  of  compressed  data.  For  HSI  data,  the  distribution  of  elements  in  the  resulting 
residual  matrix,  i.e.,  the  original  data  subtracted  by  its  low-rank  representation,  exhibits 
a  low  entropy  relative  to  the  original  data  that  favors  high-compression  ratio.  We  showed 
both  theoretically  and  empirically  that  randomized  methods  combined  with  residual-coding 
algorithms  can  lead  to  effective  lossless  compression  of  HSI  data.  We  conducted  numerical 
tests  on  real  large-scale  HSI  data  that  shows  promise  in  this  case.  In  addition,  we  showed 
that  randomized  techniques  can  be  applicable  for  encoding  on  resource-constrained  on-board 
sensor  systems,  where  the  core  matrix- vector  multiplications  can  be  easily  implemented  on 
computing  platforms  such  as  graphic  processing  units  or  field-programmable  gate  arrays. 

1.3.5  Joint  Blind  Deconvolution  and  Spectral  Unmixing  of  Hyperspectral  Im¬ 
ages  [17] 

Our  interest  here  was  spectral  imaging  for  space  object  identification  based  upon  imag¬ 
ing  using  simultaneous  measurements  at  different  wavelengths.  AMOS  sensors  can  collect 
simultaneous  images  ranging  from  visible  to  LWIR.  On  the  other  hand,  multiframe  blind 
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Figure  17:  The  reconstruction  of  compressive  hyperspectral  camera,  (a)  The  raw  measure¬ 
ment  of  the  color  candies,  (b)  The  48  spectral  channels  reconstruction  between  450nm  and 
685  nm. 

deconvolution  (MFBD)  has  demonstrated  success  by  acquiring  near-simultaneous  multiple 
images  for  reconstructing  space  objects,  and  another  success  has  been  shown  through  adding 
phase  diversity  (PD)  by  splitting  the  light  beam  in  channels  with  different  phase  functions. 
So  far,  most  MFBD  and  PD  applications  have  been  focused  on  monochromatic  images, 
with  a  few  MFBD  studies  on  multispectral  images,  also  called  the  wavelength  diversity.  In 
particular,  B.  Calef  has  shown  that  wavelength-diverse  MFBD  is  a  promising  technique  for 
combining  data  from  multiple  sensors  to  yield  a  higher-quality  reconstructed  image.  Here,  we 
presented  optimization  algorithms  to  blindly  deconvolve  observed  blurred  and  noisy  hyper¬ 
spectral  images  with  phase  diversity  at  each  wavelength  channel.  We  used  the  facts  that  at 
longer  wavelengths,  turbulence  effects  on  the  phase  are  less  severe,  while  diffraction  effects  at 
shorter  wavelengths  are  less  severe.  Moreover,  because  the  blurring  kernels  of  all  wavelength 
channels  essentially  share  the  same  optimal  path  difference  (OPD)  function,  we  have  greatly 
reduced  the  number  of  parameters  in  the  blurring  kernel.  We  modeled  the  true  hyperspectral 
object  by  a  linear  spectral  unmixing  model,  which  reduces  the  number  of  pixels  to  be  recov¬ 
ered.  Because  the  number  of  known  parameters  is  far  greater  than  the  number  of  unknowns, 
the  method  enjoys  an  enhanced  capability  of  successful  reconstruction.  We  simultaneously 
reconstructed  the  true  object,  estimated  the  blurring  kernels,  and  separate  the  object  into 
spectrally  homogeneous  segments,  each  characterized  by  its  support  and  spectral  signature. 
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Figure  18:  The  spectral  reconstruction  of  color  candies,  (a)  The  reference  image  measured  by 
a  color  camera,  (b)  to  (f)  are  the  averaged  reconstruction  spectra  of  orange,  blue,  red,  yellow 
and  green  candies,  respectively.  The  reference  spectra  are  measured  by  a  fiber  spectrometer. 

an  important  step  for  analyzing  the  material  compositions  of  space  objects. 

1.3.6  Shape  and  Pose  Recovery  of  Solar- Illuminated  Surfaces  from  Compressive 
Spectral-Polarimetric  Image  Data  [5] 

We  performed  simulation-based  studies  of  the  use  of  compressively  sensed  spectral  polari- 
metric  spatial  image  data  from  a  solar  illuminated  reflecting  surface  to  recover  its  material 
signature,  three-dimensional  (3D)  shape,  pose,  and  degree  of  surface  roughness.  The  spa¬ 
tial  variations  of  the  polarimetric  BRDF  around  glint  points  contain  unique  information 
about  the  shape  and  roughness  of  the  reflecting  surface  that  is  revealed  most  dramatically 
in  polarization-difference  maps  from  which  the  spatially  generalized  diffuse-scattering  con- 
tibutions  to  brightness  are  largely  absent.  Here  we  employed  a  specific  compressed-sensing 
protocol,  the  so-called  Coded- Aperture  Snapshot  Spectral  Polarimetric  Imager  advanced  re¬ 
cently  by  Tsai  and  Brady  under  partial  support  from  this  grant  and  discussed  in  detail 
in  Section  1.2  of  this  report,  to  simulate  noisy  measurements  from  which  these  surface  at¬ 
tributes  are  recovered  in  a  sequential  manner.  Even  in  the  presence  of  additive  sensor  noise, 
the  recovery  of  these  surface  attributes  seemed  to  be  quite  robust. 

1.3.7  Image  Reconstruction  from  Double  Random  Projection  [18] 

We  advanced  double  random  projection  methods  for  reconstruction  of  imaging  data.  The 
methods  draw  upon  recent  results  in  the  random  projection  literature,  particularly  on 
lowrank  matrix  approximations,  and  the  reconstruction  algorithm  has  only  two  simple  and 
non-iterative  steps,  while  the  reconstruction  error  is  close  to  the  error  of  the  optimal  low- 
rank  approximation  by  the  truncated  singular-value  decomposition.  We  extended  the  often- 
required  symmetric  distributions  of  entries  in  a  random-projection  matrix  to  asymmetric 
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distributions,  which  can  be  more  easily  implementable  on  imaging  devices.  Experimental 
results  were  provided  on  the  subsampling  of  natural  images  and  hyperspectral  images,  and 
on  simulated  compressible  matrices.  Comparisons  with  other  random  projection  methods 
were  also  provided. 

1.3.8  Deblurring  and  Sparse  Unmixing  of  Hyperspectral  Images  using  Multiple 
Point  Spread  Functions  [19] 

We  advanced  an  iterative  approach  to  solve  separable  nonlinear  least  squares  problems  aris¬ 
ing  in  the  estimation  of  wavelength-dependent  point  spread  function  (PSF)  parameters  for 
hyperspectral  imaging.  A  variable  projection  Gauss-Newton  method  was  used  to  solve  the 
nonlinear  least  squares  problem.  An  analysis  shows  that  the  Jacobian  can  be  potentially 
very  ill-conditioned.  To  deal  with  this  ill-conditioning,  we  used  a  combination  of  subset 
selection  and  other  regularization  techniques.  Experimental  results  related  to  hyperspectral 
PSE  parameter  identification  and  star  spectrum  reconstruction  illustrate  the  effectiveness  of 
the  resulting  numerical  scheme. 

1.3.9  Estimation  of  Atmospheric  PSF  Parameters  for  Hyperspectral  Imaging 

[20] 

A  numerical  approach  was  provided  for  deblurring  and  sparse  unmixing  of  ground-based 
hyperspectral  images  (HSI)  of  objects  taken  through  atmospheric  turbulence.  Hyperspectral 
imaging  systems  capture  a  3D  datacube  (tensor)  containing:  2D  spatial  information,  and  ID 
spectral  information  at  each  spatial  location.  Pixel  intensities  vary  with  wavelength  bands 
providing  a  spectral  trace  of  intensity  values,  and  generating  a  spatial  map  of  spectral  vari¬ 
ation  (spectral  signatures  of  materials).  The  deblurring  and  spectral  unmixing  problem  is 
quite  challenging  since  the  point  spread  function  (PSF)  depends  on  the  imaging  system  as 
well  as  the  seeing  conditions  and  is  wavelength  varying.  We  showed  how  to  efficiently  con¬ 
struct  an  optimal  Kronecker  product-based  preconditioner,  and  provide  numerical  methods 
for  estimating  the  multiple  PSFs  using  spectral  data  from  an  isolated  (guide)  star  for  joint 
deblurring  and  sparse  unmixing  the  HSI  datasets  in  order  to  spectrally  analyze  the  image 
objects.  The  methods  were  illustrated  with  numerical  experiments  on  a  commonly  used  test 
example,  a  simulated  HSI  of  the  Hubble  Space  Telescope  satellite. 
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